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Abstract 

We study the thermodynamic Casimir effect in thin films in the three 
dimensional XY universality class. To this end, we simulate the improved 
two component 0^ model on the simple cubic lattice. We use lattices up to 
the thickness Lq = 33. Based on the results of our Monte Carlo simulations 
we compute the universal finite size scaling function 9 that characterizes the 
behaviour of the thermodynamic Casimir force in the neighbourhood of the 
critical point. We confirm that leading corrections to the universal finite size 
scaling behaviour due to free boundary conditions can be expressed by an 
effective thickness -Z>o,e// = -^o + -^s, with Lg = 1.02(7). Our results are 
compared with experiments on films of ^He near the A-transition, previous 
Monte Carlo simulations of the XY model on the simple cubic lattice and field- 
theoretic results. Our result for the finite size scaling function 9 is essentially 
consistent with the experiments on films of ^He and the previous Monte Carlo 
simulations. 



Keywords: A-transition, Classical Monte Carlo simulation, thin films, thermody- 
namic Casimir effect 



1 Introduction 



In 1978 Fisher and de Gennes jlj realized that there should be a so called "thermody- 
namic" Casimir effect. This means that a force emerges when thermal fluctuations 
are restricted by a container. Thermal fluctuations extend to large scales in the 
neighbourhood of critical points. In the thermodynamic limit, in the neighbour- 
hood of the critical point, various quantities diverge following power laws. E.g. the 
correlation length, which measures the spatial extension of fluctuations, behaves as 



where t = (T — Tc)/Tc is the reduced temperature and Tc the critical temperature. 
^0,+ and ^0,- are the amplitude of the correlation length in the high and low temper- 
ature phase, respectively. While ^o,+ and ^o,- depend on the microscopic details of 
the system, the critical exponent u and the ratio ^o,+/^o,- are universal. This means 
that they assume exactly the same values for all systems within a given universality 
class. A universality class is characterized by the spatial dimension of the system, 
the range of the interaction and the symmetry of the order parameter. The modern 
theory of critical phenomena is the Renormalization Group (RG). For reviews see 
e.g. [21 [Sj |4l [5] . Here we consider the XY universality class in three dimensions with 
short range interactions. This universality class is of particular interest, since the 
A-transition of ^He is supposed to share this universality class. The most accurate 
experimental results for critical exponents and universal amplitude ratios for a three 
dimensional system have been obtained for this transition; for a review see [6]. 

The critical behaviour is modified by a confining geometry. If the system is 
finite in all directions, thermodynamic functions have to be analytic. I.e. a singular 
behaviour like eq. ([T]) is excluded. As a remnant of such singularities there remains 
a peak in the neighbourhood of the transition. With increasing linear extension 
the hight of the peak increases and the temperature of the maximum approaches 
the critical temperature. This behaviour is described by the theory of finite size 
scaling (FSS). For reviews see [Tl[S]. In general the physics in the neighbourhood of 
the transition is governed by the ratio Lq/C,, where Lq is the linear extension of the 
container and ^ the correlation length of the bulk system. Furthermore it depends 
on the geometry of the container and on the type of the boundary conditions that 
the container imposes on the order parameter. For a review on experimental studies 
of "^He near the A-transition in confining geometries see |9j. 

Here we study thin films. Thin films are finite in one direction and infinite in 
the other two directions. In this case singular behaviour is still possible. However 
the associated phase transition belongs to the two-dimensional universality class. 
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I.e. in the case of U{1) symmetry, a Kosterlitz-Thouless (KT) transition [T0l[TT| [T2] 
is expected. In [13] we have confirmed the KT-nature of this transition and have 
studied the scaling of the transition temperature with the thickness of the film. 
Recently [Hj we determined the finite size scaling behaviour of the specific heat of 
thin films. Here we investigate the thermodynamic Casimir force in thin films in 
the three dimensional XY universality class. 

From a thermodynamic point of view, the Casimir force per unit area is given 

by 

p . . - _^f£^ (2) 

^ casimir or V / 

where Lq is the thickness of the film and fex = ffUm — -^o/sd is the excess free 
energy per area of the film, where ffUm is the free energy per area of the film and 
fsD the free energy density of the thermodynamic limit of the three dimensional 
system; see e.g. [TB]. Finite size scaling predicts that the Casimir force behaves as 

k T 

Pcasimir 7^ 

where d{x) is a universal finite size scaling function. In HH [m ^He films of 
thicknesses up to 588 A have been studied. These experiments show clearly that the 
thermodynamic Casimir force is indeed present. Throughout it is negative. In the 
low temperature phase of the three dimensional bulk system it shows a pronounced 
minimum. The data are essentially consistent with the prediction eq. ([3]). The 
minimum of 0{x) is located at a; = t{Lo/^oy^'^ ^ —5.5. 

It has been a challenge for theorists to compute the finite size scaling func- 
tion 6{x). Krech and Dietrich [TU] have computed it in the high temperature 
phase using the e-expansion up to 0(e). This result is indeed consistent with the 
measurements on ^He films. Deep in the low temperature phase, the spin wave 
approximation should provide an exact result. It predicts a negative non- vanishing 
value for 6{x). However the experiments suggest a much larger absolute value for 
6{x) in this region. Until recently a reliable theoretical prediction for the minimum 
of 6{x) and its neighbourhood was missing. Using a renormalized mean-field ap- 
proach the authors of [20] have computed 6{x) for the whole temperature range. 
Qualitatively they reproduce the features of the experimental result. However the 
position of the minimum is by almost a factor of 2 different from the experimental 
one. The value at the minimum is wrongly estimated by a factor of about 5. 

^Following the literature, in cq. fS]) wc shall use ^o.+ in the following. 
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Only quite recently Monte Carlo simulations of the XY model on the simple 
cubic lattice [211 1221 [23] provided results for 9{x) which essentially reproduce the 
experiments on ^He films [T6l [T7] . These simulations were performed with lattices 
of a thickness up to Lq = 16 [22] and up to Lq = 20 [23]. The authors of [23] 
pointed out that for these lattice sizes corrections to scaling still play an important 
role. The purpose of the present work is to get accurate control over the leading 
corrections to scaling, allowing us to compute 6{x) with reliable error bars. 

As the first step in this direction we simulate the improved two-component 
0^ model instead of the XY model. For a precise definition of these models see 
the section below. This way we avoid leading corrections to scaling which are 
oc Lq^ with uj = 0.785(20) [23]; similar values for the exponent are obtained with 
field-theoretic methods; for a review see e.g. [5j. In order to mimic the vanishing 
order parameter that is observed at the boundaries of "^He films, Dirichlet boundary 
conditions with vanishing field are imposed. These lead to corrections oc Lq ^ [25j. 
These corrections can be eliminated by replacing the thickness Lq by an effective 
one Loeff = Lq + Ls, where Lg = 1.02(7) [13] for the model that we have simulated. 
I 

This paper is organized as follows: First we define the model and the observables 
that we have measured. Next we discuss the finite size scaling behaviour of the 
Casimir force. In particular, we discuss corrections to scaling caused by the Dirichlet 
boundary conditions. We outline the method used to compute the Casimir force. 
We discuss the simulations that have been performed and analyze our data. We 
compare our results with experiments [161 [U]) previous Monte Carlo simulations 
[22l [23] and the e-expansion [TH [19]. Finally we summarize and conclude. 

2 The model and the observables 

We study the two component 0^ model on the simple cubic lattice. We label the 
sites of the lattice by ). The components of x might assume the 

values Xi G {1, 2, . . . , Lj}. We simulate lattices of the size Li = L2 = L and 
Lq <^ L. In 1 and 2-direction we employ periodic boundary conditions and free 
boundary conditions in 0-direction. This means that the sites with xq = 1 and 
Xq = Lq have only five nearest neighbours. This type of boundary conditions could 
be interpreted as Dirichlet boundary conditions with as value of the field at = 

^ In the literature, replacing Lq by L^^eff ^ Lq + Lg to account for surface corrections, was 
first discussed by Capehart and Fisher [26] in the context of the surface susceptibility of Ising 
films. 
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and xq = Lq + I. Note that viewed this way, the thickness of the film is Lq + 1 rather 
than Lq. This provides a natural explanation of the result Lg = 1.02(7) obtained in 
[T3] and might be a good starting point for a field theoretic calculation of Lg. The 
Hamiltonian of the two component (p^ model, for a vanishing external field, is given 
by 



<x,y> 



where the field variable (px is a vector with two real components. < x,y > denotes 
a pair of nearest neighbour sites on the lattice. The partition function is given by 



n 



1) 



exp{-n). (5) 



Note that following the conventions of our previous work, e.g. [28] , we have absorbed 
the inverse temperature /3 into the Hamiltonian. [f] In the limit A — > oo the field 
variables are fixed to unit length; i.e. the XY model is recovered. For A = we get 
the exactly solvable Gaussian model. For < A < cxd the model undergoes a second 
order phase transition that belongs to the XY universality class. Numerically, 
using Monte Carlo simulations and high-temperature series expansions, it has been 
shown that there is a value A* > 0, where leading corrections to scaling vanish. 
Numerical estimates of A* given in the literature are A* = 2.10(6) [27J, A* = 2.07(5) 
[28] and most recently A* = 2.15(5) [21]. The inverse of the critical temperature 
f3c has been determined accurately for several values of A using finite size scaling 
(FSS) [23]. We shall perform our simulations at A = 2.1, since for this value of A 
comprehensive Monte Carlo studies of the three-dimensional system in the low and 
the high temperature phase have been performed [TSj [231 1291 [3Q|. At A = 2.1 one 
gets Pc = 0.5091503(6) [23]- Since A = 2.1 is not exactly equal to A*, there are 
still corrections oc L~^, although with a small amplitude. In fact, following [23], it 
should be by at least a factor 20 smaller than for the standard XY model. 

In [13] we find for A = 2.1 by fitting the data for the second moment correlation 
length in the high temperature phase 

^2nd = 0.26362(8)^-°-^^^^ X [1 + 0.039(8)t°-^^^ - 0.72(4)t] , (6) 

where t = 0.5091503 — /5. Here we shall use u = 0.6717(1) [23] as value of the critical 
exponent of the correlation length. Recent experiments on the A-transition of "^He 



^Therefore, following [4] we actually should call it reduced Hamiltonian. 
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suggest a slightly smaller value: v = 0.6709(1) [3T]. This discrepancy is however 
not crucial for the present study. Note that in the high temperature phase there 
is little difference between C,2nd and the exponential correlation length ^exp which is 
defined by the asymptotic decay of the two-point correlation function. Following 

m- 

lim^ = 1.000204(3) , (t > 0) (7) 

C,2nd 

for the thermodynamic limit of the three-dimensional system. 

2.1 The internal energy and the free energy 

The reduced free energy density is defined as 

f = -Y^\ogZ. (8) 

I.e. compared with the free energy density /, a factor fc^T is skipped. 

Note that in eq. (jl]) j3 does not multiply the second term. Therefore, strictly 
speaking, (3 is not the inverse of fc^T. In order to study universal quantities it is 
not crucial how the transition line in the f3-X plane is crossed, as long as this path is 
not tangent to the transition line. Therefore, following computational convenience, 
we vary f3 at fixed A. Correspondingly we define the (internal) energy density as 
the derivative of the reduced free energy density with respect to (3. Furthermore, 
to be consistent with our previous work [H], we multiply by —1: 

(9) 

It follows 




which can be easily determined in Monte Carlo simulations. From eqs. (p[9i) it 
follows that the free energy density can be computed as 

/(/?) = /(/5o) - l^dpEiP) . (11) 
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3 The finite size scaling behaviour of the thermal 
Casimir force 



Let us discuss the scaling behaviour of the reduced excess free energy. Since we 
study an improved model we ignore corrections cx Lq in the following. We take 
into account leading corrections due to the boundary conditions by replacing the 
thickness Lq of the film by Lo,e// = Lo + Ls at the appropriate places. We split the 
free energies in singular (s) and non-singular (ns) parts: 

fex(t,Lo) = ffilmit, Lq) — Lofsoit) 

= ffUm,s(t, Lq) + LQ f,ff^nsfns{i) ~ Lof^D^t) — Lo/nsit) 

= L-l^fh{x)+L,hD,s{t)+Lnsfns{t) (12) 

where h{x) = Ll^^j.j[ffiim,s{t,Lo) - Lo,e///3D,s(i)] is a universal finite size scaling 
function and x = t[Lo^eff/^oY^'^- Following RG theory the non-singular part is not 
affected by finite size effects. However it is not clear a priori how Dirichlet boundary 
conditions affect the non-singular part of the free energy. Therefore we allow for 
Lns = LQ eff,ns — Lq and Lns 7^ Lg. Taking the derivative with respect to Lq we 
get the thermodynamic Casimir force per area [15] 

PFcasimiT = ^ = 2Lq gjy/i(x) — L^^^^ ^—xh\x) = Lq^j-j6{x) (13) 

where 6{x) = 2h{x) — -xh'{x). 



4 Computing the Casimir force on the lattice 

Here we follow essentially the approach of [22]. For an alternative method see 
[2T| [23] . On the lattice the thickness Lq assumes integer values. Therefore we 
approximate the derivative for half-integer values of Lq as 



^ Af,M Lq) = fi/3, Lq + 1/2) - fif3, Lq - 1/2) + f,Di/3) ■ (14) 

L=Lo 



dL 

Correspondingly we define 

AE,,{(3, Lq) = E{(3, Lq + 1/2) - E{(3, Lq - 1/2) - E^DiP) (15) 
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where E{f], Lq) is the energy per area of a thin film and E^oiP) the energy density 
of the three dimensional system. Analogous to eq. (ITT]) we compute 

AU{(3) = - / d(3AE,,{P) (16) 

where /3o is chosen such that ^{Po) ^ Lq and hence the Casimir force vanishes. 

5 Numerical results 

In [13] we have studied the specific heat of thin films in the two component (f)^ 
model at A = 2.1. To this end, we had determined the energy density for the three 
dimensional thermodynamic limit and for films of the thicknesses Lq = 8, 16 and 
32 for a large number of /3-values. In order to compute the derivative with respect 
to Lq, we have complemented these simulations by ones for Lq = 9, 17 and 33. The 
Monte Carlo algorithm that has been used is the same as in [T^: An update-cycle 
is composed of a Metropolis sweep, a few overrelaxation sweeps and single ^32] and 
wall ^33j cluster updates. One sweep means that a local update is performed at 
each site of the lattice ones. As random number generator we have used the SIMD- 
oriented Fast Mersenne Twister algorithm [33] • In table [1] we have summarized the 
statistics of our runs. In total these simulations took about 3 years of CPU-time 
on a single core of a Quad-Core Opteron(tm) 2378 CPU (2.4 GHz). 

Using these data, we have computed AEe^ for Lq = 8.5, 16.5 and 32.5. One 
should note that the statistical error of E{P,Lq + 1/2) — E{P,Lq — 1/2) is much 
larger than that of E3d{(3). In order to obtain Afe^ we have numerically integrated 
AEex using the trapezoidal rule: 



where /3j are the values of (3 we have simulated at. They are ordered such that 
Pi+i > Pi for all i. We have chosen Pq = 0.49, 0.5 and 0.505 for Lq = 8.5, 16.5 and 
32.5. We find that AE^x is equal to zero within error bars up to values of P that 
are slightly larger than the Pq that we have chosen. 

The estimate obtained from the integration is affected by statistical and sys- 
tematical errors. The statistical one can be easily computed, since the AE^^ are 
obtained from independent simulations: 



AfexiPn) ~ 2 ^^^+1 ~ {AE,x{P^+l) + AE,x{P^)) 



(17) 



j=0 



e\-Af,x{Pn)) 



{Pl-Po? 



e^[AEUP^)] + 




e^[AE,x{Pn)] 



4 
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Table 1: We characterize our new simulations. In the first column we give the 
thickness Lq of the film. In the second column we give the linear size L = Li = L2 
of the lattice in the other two directions. In the third and fourth column we give 
the upper and lower bound of the interval in j3 that has been simulated. In the 
fifth column we give the step size A/5 that we used. E.g. Prnin = 0.49, Pmax = 0.519 
and Ap = 0.001 means that (3 = 0.49, 0.491, 0.492, . . ., 0.519 have been simulated. 
Finally, in the last column we give the number of measurements (stat) that we have 
performed for each of the simulations. 





Li — L2 




Pmax 


AP 


stat 


9 


64 


0.49 


0.519 


0.001 


5 X 10^ 


9 


128 


0.52 


0.527 


0.001 


2 X 10^ 


9 


256 


0.528 


0.56 


0.001 


10^ 


9 


256 


0.562 


0.58 


0.002 


10^ 


9 


512 


0.536 


0.539 


0.001 


10^ 


9 


512 


0.5395 


0.548 


0.0005 


10^ 


9 


512 


0.548 


0.57 


0.001 


10^ 


9 


1024 


0.539 


0.548 


0.0005 


10^ 


17 


256 


0.527 


0.55 


0.001 


2 X 10^ 


17 


512 


0.5 


0.512 


0.001 


10^ 


17 


512 


0.5125 


0.529 


0.0005 


10^ 


17 


512 


0.53 


0.55 


0.001 


10^ 


17 


1024 


0.5205 


0.529 


0.0005 


8 X 10^ 


33 


256 


0.502 


0.50875 


0.00025 


4 X 10^ 


33 


512 


0.509 


0.5128 


0.0002 


3 X 10^ 


33 


1024 


0.513 






10^ 


33 


1024 


0.5132 






8 X 10^ 
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+ 4 ^'[^E^M)] (18) 

i=l 

where denotes the square of the statistical error. 

In order to estimate the error due to the finite step size — we have redone 
the integration, skipping every second value of /?; i.e. doubling the step size. In 
all three cases (i.e. Lq = 8.5, 16.5 and 32.5), the results were consistent within the 
statistical errors. Therefore we are confident that the systematical error due to the 
finite step size is smaller than the statistical one. 

In figure[I]we have plotted —L^Afex as a function of —t[Lo/^oY^'^ , where we have 
used p = 0.6717 and = 0.26362, eq. (Ej). We find that throughout the function 
assumes a negative value. In all cases it has a single minimum at t[Lo/^o]^'''^ ~ 
—5. The position of the minimum j3min{Lo) can be easily determined: It is given 
by the zero of AE^x- We have computed PminiLo) by linearly fitting AE^x in 
the neighbourhood of the minimum. In addition to Lq = 8.5, 16.5 and 32.5 we 
performed simulations for Lq = 6.5, 7.5, 9.5, 12.5 and 24.5 at a few values of /3 in 
the neighbourhood of Pmin- Our results are summarized in table [2l 

The curves for Lq = 8.5, 16.5 and 32.5 plotted in figured] do not fall on top of 
each other. E.g. both the position and the value of the minimum are quite different 
for different Lq. In order to take corrections into account we have replaced Lq by 
Lo,eff = Lq + Lg, where Lg = 1.02(7) [13]. To this end, in figure [2] we have plotted 
— Lp g^jA/e^,. as a function of — t[Lo,e///Co]^''^! where we have used the central value 
of the shift Lg = 1.02. Now indeed the distance between the curves for different Lq 
is much reduced compared with figure [H The results for Lq = 16.5 and Lq = 32.5 
are almost consistent within error bars. Note that using Lg = 0.95 the matching of 
the data for different Lq seems to be better than for Lg = 1.02. 

Let us discuss in more detail the results obtained for the minimum of the finite 
size scaling function 6{x). Using the numbers given in the third column of table [2] 
and Lg = 1.02 we get -Afex,minLl^ff = -1.365(3), -1.341(6) and -1.311(19) for 
Lq = 8.5, 16.5 and 32.5, respectively. Using instead Lg = 0.95 we get —1.335(3), 
— 1.325(6) and —1.302(19) for Lq = 8.5, 16.5 and 32.5, respectively. As our final 
result we take the one obtained from Lq = 32.5 and Lg = 1.02: 

= -1.31(3) , (19) 

where the error that is quoted takes into account the statistical error and the un- 
certainty oi Lg. 

Next we have fitted our results for jSmm with the ansatz 

(20) 
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Figure 1: We plot —Li^^fex as a function of —t^LQ/^Qy/" for Lq = 8.5, 16.5 and 
32.5, where we use u = 0.6717 and = 0.26362. For a discussion see the text. 
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Figure 2: We plot —Ll^^jj^Af^x as a function of — t(Lo,e///^o)^^'^ for -^o = 8-5, 
16.5 and 32.5, where we use u = 0.6717, = 0.26362 and Lo^^ff = Lq + with 
Lg = 1.02. For a discussion see the text. 
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Table 2: We give the position (3min of the minimum of the Casimir force and its 
value — A/ea;. function of the thickness Lq. 







ex, min 


6.5 


0.54432(2) 




7.5 


0.53814(2) 




8.5 


0.53354(2) 


-0.001582(3) 


9.5 


0.53010(2) 




12.5 


0.52348(2) 




16.5 


0.51886(2) 


-0.0002494(11) 


24.5 


0.51463(2) 




32.5 


0.51279(2) 


-0.0000348(5) 



Table 3: We fit our results for Pmin with ansatz ( f20|) . Lo,mm is the smallest thickness 
of the film that is included into the fit. 



Lo^min 


•^min 


c 


xVd.o.f. 


6.5 


-4.942(6) 


1.20(4) 


1.47 


7.5 


-4.945(8) 


1.17(7) 


1.72 


8.5 


-4.956(10) 


1.04(10) 


1.32 


9.5 


-4.952(12) 


1.10(12) 


1.64 



where t^in = - Pmin- We have used u = 0.6717, = 0.26362, (3^ = 0.5091503 
and Lg = 1.02 as input and c and Xmin as parameters of the fit. The term (1 + 
c-tmin) parametrizes analytic corrections. We did not include a correction with 
the exponent 6' ~ 1.2 [3S|, since within the accuracy of our data they can not 
be discriminated from the leading analytic correction. The results of these fits are 
summarized in tableO The x^/d.o.f. is reasonably small starting from Lo,mm = 6.5, 
where all thicknesses Lq with Lq > LQ^min are included into the fit. Also the 
estimates for Xmin and c do not change much as Lo,mm is changed. In order to 
check the dependence of the results on the value of Lg we have repeated the fit for 
LQ,min = 8.5 using Lg = 0.95 instead of Lq = 1.02. We get the Xmin = —4.943(10), 
c = 0.70(10) with x^/d.o.f.= 1.12. As our final result we take 

x^in = -4.95(3) (21) 
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where the error bar covers the statistical error and the uncertainty of Lg. 

Our result for x^m can be compared with those given in the literature. The 
experimental works [16] give Xmin = —5.45(12) (no final result for 9min is quoted) 
and [I7j Xmin = —5.7(5) and 9min = —1.30(3). In pjj and [IT] the convention 
X = tL\l^ is used. In order to convert to x = t[Lo/^o]^^'' we have used .^o = 1-422 A 
for ^He at vapour pressure as discussed in section 4.2 of [Ti] . 

In his Monte Carlo study of the XY model on the simple cubic lattice, Hucht 
[22] finds x = —5.3(1) and 9min = —1.35(3). The authors of [23| have used two 
different ansatze for the corrections to scaling. Using the first one, they arrive at 
Xmin = —5.43(2) and 9min = —1.396(6) and using the second one at Xmin = —5.43(2) 
and 9min = —1.260(5). Note that the ansatze used by [23], see eqs. (20,21,22,23) of 
[23], provide only an overall x- independent factor; therefore they do not allow for 
any correction to scaling of Xmin- 

Our results for 9min is in good agreement with both the experiment \V7l as well 
as with previous Monte Carlo studies of the XY model [22l[23]. On the other hand, 
the position of the minimum Xmin differs by several times the quoted error bar from 
both the experiment [Ml Il7| as well as from Monte Carlo studies of the XY model 
[221 [23]. 

Finally, in figure [3] we take into account the analytic corrections that we have 
detected fitting the position tmin of the minimum of the Casimir force. To this 
end, we have replaced the argument t{LQ^eff/ioY^'' by t{l + l.{iAt){LQ^ef f / ioY^" ■ 
The coefficient of the analytic correction is taken from the fit where we have fixed 
Ls = 1.02 and L^min = 8.5. Now we find an almost perfect match between the 
curves obtained from Lq = 8.5, 16.5 and 32.5. 

5.1 Comparison with other theoretical approaches 

Krech and Dietrich [TH [19] have computed the finite size scaling function 9 in the 
high temperature phase using the e-expansion up to 0(e). In figure [H we plot their 
result for the XY universality class (A^ = 2) setting e = 1. For comparison we plot 
our results for Lq = 8.5, Lq = 16.5 and Lq = 32.5. We have taken into account 
leading boundary corrections by replacing Lq by L^ eff = Lq + L^, where we have 
taken Lg = 1.02. 

Comparing with the e-expansion we can estimate the systematical error caused 
by setting A/ex (/So) = in eq. ( ITTll : We read off from the e-expansion that 9 ~ 
-0.0035, -0.0022, -0.0015 for our choices of po for Lq = 8.5, 16.5 and 32.5. Taking 
into account this error, we see a good agreement of our Monte Carlo results and the 
e-expansion down to Lq/C, ~ 1. The curve obtained from the e-expansion flattens 
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Figure 3: We plot — Lq ^j-j^Afex is plotted as a function of -t{l + lMt){Lo^eff /^oY^" 
for Lq = 8.5, 16.5 and 32.5, where we use z/ = 0.6717, = 0.26362 and Lo^^ff = 
Lq + Ls with Ls = 1.02. For a discussion see the text. 
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as the critical point is approached. At the critical point the slope vanishes. In 
contrast, in our case, the curve steepens as the critical point is approached. 

The authors of [20] have computed the finite size scaling function 9 using a 
renormalized mean field approach. While this approach correctly reproduces qual- 
itative features of the finite size scaling function 9 it fails to give quantitatively 
accurate results. In particular the authors of |20j find Xmin = — tt^ ^ —9.8696 and 
^mm ~ —6.92. I.e. The position of the minimum is overestimated by about a factor 
of 2 and its value by a factor of about 5. 

The authors of [36] incorporate fiuctuation effects into their mean-field analysis. 
Qualitatively they get the finite size scaling function 9 for the whole range of x 
right. They adjust the two parameters of their solution for the low temperature 
phase (eq. (17) of [36] ) such that the minimum of 9 found in the experiments 
[T6t [T7] is reproduced. We did the same exercise, adjusting to our result for the 
minimum of 9. We find no good match between 9 computed in ^36j and ours. Our 
minimum is much more peaked than that of 



5.2 Comparison with experimental results 

Finally we compare our result for the finite size scaling function 9 with experiments 
[T6t [T7] . In [16] films of thicknesses between 298 and 588 A have been studied. In 
figure [5] we have plotted the data obtained from capacitor 1 which corresponds to 
the thickness 575A of the film at temperatures T > T\. This set of data is the 
smoothest among the five sets given in [161 EZ]- The results of [17] are, in the 
range of temperatures we are interested in, less precise than those of [I6]. In the 
tables provided by the authors [37] the finite size scaling function 9 is given as a 
function of {T/Tx — . In order to compare with our results we have converted 

this to (T/Tx — l)(-^o/^o)^^'') using = 1.422A. For comparison we give our result 
obtained from Lq = 16.5, where we have taken into account the boundary correction 
by replacing Lq by Lo^eff and the leading analytic correction as discussed above. 
Furthermore, we give the asymptotic value [381 



lim 9{x) = ^ -0.04783 (22) 

obtained from the spin wave approximation. 

As already observed by the authors of [22l [23] there is a qualitative agreement 
among the result obtained from Monte Carlo simulations of lattice models and the 
experiment. There is a reasonable agreement of the position of the minimum Xmin, 
as already discussed above. For x < Xmin our result is in good agreement with 
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Figure 4: We plot the finite size scaling function as a function of Lq/^ in the high 
temperature phase obtained by Krech and Dietrich fiSl [19] using the e-expansion. 
For comparison we give our results obtained for Lq = 8.5, Lq = 16.5 and Lq = 32.5. 
In the case of our data the leading boundary correction is taken into account by 
replacing Lq by i^o,e// = Lq + Ls with Lg = 1.02. For a discussion see the text. 
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that of the experiment. In contrast for x > Xmin the experimental value of 6 is 
clearly smaller than ours. For x ~ —30 the experimental result (capacitor 1 of [T6] ) 
assumes ~ —0.19 and is decreasing again for smaller values of x. This is clearly 
different from the prediction ( |22l) of the spin wave approximation. 

In ref. [39] it was argued that this discrepancy could be explained by fluctuations 
of the surface resulting in 

lim e{x) = ^ -0.1315 . (23) 

This goes indeed in the right direction, can however not fully explain the difference 
between the experimental result and the theoretical prediction fl22]) . 

6 Summary and Conclusions 

We have simulated the improved two component (f)'^ model on the simple cubic 
lattice. This model shares with the A-transition of ^He the three dimensional XY 
universality class. We consider the thin film geometry. In order to mimic the 
vanishing order parameter at the surface of ^He films near the A-transition, we 
impose Dirichlet boundary conditions with a vanishing field. 

Restricting the system to a finite geometry leads to an effective force called 
thermodynamic Casimir force. Its behaviour in the neighbourhood of the critical 
point is characterized by a universal finite size scaling function. We have computed 
this function and have compared our result with that obtained from experiments 
on films of ^He [T6l [T7l [37] , previous Monte Carlo simulations of the XY model on 
the simple cubic lattice [2T], [22l [23] , field theoretic methods [H], [19] and mean field 
approaches [201 [3S] • 

The thermodynamic Casimir force is given as minus the derivative of the excess 
free energy of the film with respect to its thickness Lq. On the lattice, this is 
approximated by the finite difference of films of the thickness Lq + 1 /2 and Lq — 1 /2 , 
where Lq + 1/2 is integer. 

In general it is impossible to compute free energies from a single Monte Carlo 
simulation. To circumvent this problem, divide and conquer strategies are em- 
ployed. Different strategies have been proposed in [211 123] and [22] . Here we 
essentially follow [22]: We compute the derivative of the excess energy with respect 
to Lq for a dense grid of temperature values in the neighbourhood of the critical 
point. The corresponding result for the free energy is then obtained by numerical 
integration. 
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Figure 5: We plot the finite size scaling function 6{x) obtained from an experiment 
on a thin film of ^He ^ [37], where x = {T/Tx - l)(Lo/^o)^/'' with = 1.422A. 
For comparison we give our result obtained from Lq = 16.5, where we have taken 
into account boundary and analytic corrections as in figure [3l Furthermore we give 
the asymptotic value —0.04783 for x — — oo (dashed line). The authors of [20] 
have argued that fluctuations of the surface of the ^He film gives an additional 
contribution to the Casimir force. The corresponding asymptotic value —0.1315 is 
given by the dotted line. For a discussion see the text. 
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As discussed in ref. [23], corrections to scaling are numerically quite large for 
the thicknesses that can be studied at present. The main purpose of the present 
work is to get better control over these corrections as in previous work [22l [23] . 

To this end, we have studied the 0^ model at A = 2.1. In this model leading 
corrections to scaling (and finite size scaling) which are oc L^"^, with uj = 0.785(20), 
are suppressed at least by a factor of 20 compared with the XY model [21]. 

Boundary effects lead to corrections oc L^^. These corrections can be cast into 
the form Lo^ff = Lq + Lg. In [13] we have determined Lg = 1.02(7) for the two 
component 0^ model at A = 2.1 by using a finite size scaling study at the critical 
point of the three dimensional system. We have verified that this choice of Ls 
indeed eliminates the leading boundary correction in the scaling of the temperature 
of the Kosterlitz-Thouless transition [13] and the specific heat of thin films [13]. 
Also here we confirm that corrections can be essentially eliminated by replacing Lq 
by i^o,e// = Lq + Lg with Lg = 1.02(7). Remaining discrepancies can be fitted by 
analytic corrections. 

Essentially we confirm the results obtained by previous Monte Carlo simulations 
of the three dimensional XY model [22l |23] for the finite size scaling function 6. 
The main discrepancy with these previous works is the position of the minimum 
Xmin = —4.95(3) compared with Xmin = —5.3(1) [22] and Xmin = —5.43(2) [23] . 

We should note that the Monte Carlo simulation of lattice models is at the 
moment the only theoretical method that allows for a quantitatively accurate cal- 
culation of 6 in the low temperature phase not too far from the critical point. The 
e-expansion gives correctly the behaviour in the high temperature phase. The spin 
wave approximation gives the exact result in the limit x ^ —oo. The mean field 
calculation of [20] reproduces only qualitatively the features of the scaling func- 
tion. Quantitatively it is not satisfactory: the position of the minimum is wrongly 
estimated by a factor of almost 2 and its value by a factor of about 5. 

Qualitatively the Monte Carlo studies of lattice models nicely reproduce the 
finite size scaling function obtain from the experimental data [TH |T7l |37| for films 
of ''He. For x > Xmin there is a very good quantitative agreement between the 
two. In contrast, for x < Xmin the value obtained from the experiment is clearly 
smaller than that of the Monte Carlo studies. At large values of —x the spin wave 
approximation should become exact. Also in this regime, the experiments produce 
numbers that are too small compared with the theoretical one. The authors of [20] 
explain this discrepancy by fluctuations of the surface of the ^He film. Their result 
indeed reduces but not completely eliminates the difference between experiment 
and theory. Therefore it might be interesting to perform Monte Carlo simulations 
of a lattice model that incorporates fluctuations of the surface of the film. 
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